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Abstract 

We present a kinetic Monte Carlo study of the dynamical response of a Ziff-Gulari-Barshad 
model for CO oxidation with CO desorption to periodic variation of the CO pressure. We use a 
square- wave periodic pressure variation with parameters that can be tuned to enhance the catalytic 
activity. We produce evidence that, below a critical value of the desorption rate, the driven system 
undergoes a dynamic phase transition between a CO2 productive phase and a nonproductive one 
at a critical value of the period and waveform of the pressure oscillation. At the dynamic phase 
transition the period-averaged CO2 production rate is significantly increased and can be used as 
a dynamic order parameter. Wc perform a finite-size scaling analysis that indicates the existence 
of power-law singularities for the order parameter and its fluctuations, yielding estimated critical 
exponent ratios ^/u 0.12 and 7/1/ 1.77. These exponent ratios, together with theoretical 
symmetry arguments and numerical data for the fourth-order cumulant associated with the transi- 
tion, give reasonable support for the hypothesis that the observed nonequilibrium dynamic phase 
transition is in the same universality class as the two-dimensional equilibrium Ising model. 

PACS numbers: 64.60.Ht, 82.65.+r, 82.20.Wt, 05.40.-a 
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I. INTRODUCTION 



The study of nonequilibrium statistical models is a subject of current interest in a broad 
range of fields such as chemical reactions, fluid turbulence, chaos, biological populations, 
growth-deposition processes, and even economics Ql In particular, the study of surface 
reaction systems has received a great deal of attention [3|. These systems not only constitute 
a fruitful laboratory for exploring critical phenomena associated with out-of-equilibrium 
statistical physics, but they can also play an important role in the development of more 
efficient catalytic processes. Catalytic reactions are widespread in nature and have many 
industrial and technological applications 

Several experiments show that it is possible to increase the efficiency of catalytic reactions 
by subjecting the system to periodic external forcing j^O, 3,0, ^- Monte Carlo simulations 



of the Ziff, Gulari, and Barshad (ZGB) model without desorption jl^ indicate that an 
enhancement of the catalytic activity is observed when the system is perturbed by a periodic 
force that drives it briefly into the CO poisoned state However, it is well known that 

the ZGB model does not reproduce several important aspects of catalytic processes, such as 
lateral diffusion Q| and CO desorption [l^. In the present study we neglect diffusion and 
concentrate on the effects of CO desorption. For brevity we will refer to this model as the 
ZGB-k model 

Other work by the present authors jlT*] indicates that near the coexistence line between 
the active and the CO poisoned regime, the decay times of the metastable states are different 
if the ZGB-k model is driven into the CO poisoned state from the active phase, or if it is 
driven into the active phase from the CO poisoned state. Based on this result, we expect 
that the catalytic activity of the system will increase when the system is subjected to peri- 
odic variation of the external CO pressure, only when one takes into account that the time it 
takes the system to decontaminate is different from the time it takes it to contaminate. Fur- 
thermore, we show that the ZGB-k model, driven by an oscillating CO pressure , un dergoes 



a dynamic phase transition , si milar to the one observed in Ising 



anisotropic Heisenberg 
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3 lis, 



2^m 
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28| . and XY models [29|, driven by an oscillating applied fleld. 
The rest of this paper is organized as follows. In Sec. II we deflne the model and describe 
the Monte Carlo simulation techniques used. In Sec. Ill we present and discuss the numerical 



results obtained when subjecting the model to a periodic variation of the external CO 
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pressure. Finally, we present our conclusions in Sec. IV. 



II. MODEL AND SIMULATIONS 

The original ZGB model without desorption Q| describes some kinetic aspects of the 
reaction CO+0 CO2 on a catalytic surface in terms of a single parameter: the probability 
y that the next molecule arriving at the surface is CO. This parameter is proportional to the 
partial pressure of CO and will loosely be referred to as the CO pressure. The model exhibits 
two transitions, a continuous one at low CO pressure to an oxygen-saturated surface, and 
a discontinuous one at a higher CO pressure to a CO saturated surface. When desorption 
is included, the distinction between the high and low CO coverage phases disappears for a 
desorption rate k above a critical value, kc ~ 0.0406 jiol . 0]. There is evidence that the 
nonequilibrium phase transition at kr is in the same universality class as the two-dimensional 
kinetic Ising model at equilibrium . Below kc there are well-defined low and high coverage 
phases, separated by a discontinuous, nonequilibrium phase transition. This picture is con- 
sistent with experimental observations on catalytic surfaces that show transitions between 
low and high coverage phases 0, [l^ 11. 1^. 

The ZGB model with CO desorption is simulated on a square lattice of side L that 
represents the catalytic surface. A Monte Carlo simulation generates a sequence of trials: 
adsorption with probability 1 — A; and desorption with probability k. A site is selected at 
random. In desorption, if the site is occupied by a CO it is vacated, if not the trial ends. In 
adsorption, if the site is occupied the trial ends, if not a CO or O2 molecule is selected with 
probability y or 1 — y (the relative impingement rates of CO and O2), respectively. A CO 
molecule can be adsorbed at the empty site if none of its nearest neighbors are occupied by 
an O atom. Otherwise, one of the occupied neighbors is selected at random and removed 
from the surface, liberating a CO2 molecule. O2 molecules require a nearest-neighbor pair 
of vacant sites to adsorb. Once an O2 molecule is adsorbed, it dissociates into two O atoms. 
If an O atom is located next to a site filled with a CO molecule, they react to form a CO2 
molecule that escapes, leaving two sites vacant. This process mimics the CO + O ^ CO2 
surface reaction. 
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FIG. 1: Average values of the CO and O coverages, ^co a-nd Oq, and of the CO2 production rate, 
i?C02 5 shown as functions of the stationary apphed CO pressure y, for L = 100 with k = 0.02. A 
continuous, nonequihbrium phase transition occurs at yi, and a discontinuous one at 2/2 (^)- 

III. RESULTS 

Our simulations were performed on a square lattice of LxL sites, assuming periodic 
boundary conditions. The time unit is one Monte Carlo step per site (MCSS), in which each 
site is visited once, on average. 

The coverages 6*00 and 9o are defined as the fraction of surface sites occupied by CO 
and O, respectively, and Rco2 is defined as the rate of production of CO2. In Fig. [T] we 
show the dependence on the external constant CO pressure y of the average value of the 
coverages and the CO2 production rate: (6*00)5 (6*0) and (-RCO2)) respectively. There are 
two inactive regions, y < Hi {yi seems to be fairly independent of k, as expected because of 
the absence of CO on the surface at this point) and y > y2{k), corresponding to the cases in 
which the surface is saturated with O and CO, respectively. For yi < y < 2/2(^)5 there is an 
active window where the system produces CO2. Notice that the maximum value of {RCO2) 
is reached as 2/2 (^) is approached from the low-6'co phase. 

One way to perturb the catalytic oxidation process in order to increase the CO2 produc- 
tion is by switching the external pressure back and forth across the discontinuous transition 
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y2{k) . It has been shown that when the ZGB system without desorption is 

quickly oscillated between the low and high 6*00 phases by a periodic variation of the exter- 
nal pressure, there is considerable enhancement of the catalytic activity The period 
and amplitude of the driving pressure must be calibrated very carefully in order to avoid 
driving the system irreversibly toward the CO poisoned state. 

In other work we have calculated the lifetimes associated with the decay of the 
metastable states of the ZGB-k model. The system was prepared in the low (high) CO 
coverage phase with an initial pressure i/i < 1/2 (^) {Ui > y2{k)), and then y was suddenly 
changed to ?// > y2{k) {yf < y2{k)). We then measured the time it took the system to leave 
the metastable state in both cases. We found that the lifetimes depend on the direction of 
the process, the decontamination time (from high to low CO coverage) being different 
from the poisoning time Tp (from low to high CO coverage). Inspired by this result, we 
decided to subject the system to an oscillating pressure y{t) that in a period T = td + tp 
takes the values , 

{yi during the time interval td 
(1) 
yh during the time interval tp , 

located at both sides of the transition point, i.e., yi < y2{k) < yh (see Fig. |2fa)). We found 
that, for each selection of yi and yu, by tuning td and tp, the times that the driving force 
spends in the low and high coverage regions respectively, we could increase the productivity 
of the system. In Fig I21b) it is seen that the response to the periodic pressure shown in 
FiglSfa), the CO2 production rate Rco2i ^-Iso exhibits an oscillatory behavior. We therefore 
use its period- averaged value, defined as 

r = ^ ^ Rco,it)dt , (2) 

as the dynamic order parameter. For the parameters used in Fig. |21 the long-time average of 
r is (r) = 0.2683, 11% higher than the maximum average CO2 production rate for constant 
y, (i?co2)max = 0.2414. (Compare Fig. Hand FigEfb).) 

By averaging r over many periods of oscillation (of the order of 2 x 10^), we found that, 
depending on the values of td and tp, the system has two well-defined regimes: a productive 
one with (r) > 0, and a nonproductive one with (r) ^0, separated by a transition line in 
the (tp, td) plane. Fig. El shows density plots of (r) in terms of td and tp for two choices of yi 
and yh- 
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FIG. 2: (a) Applied periodic pressure of CO, y{t), that takes the values yi = 0.5 and = 0.535 
during the time intervals = 10 and tp = 150, respectively, (b) Response of the production rate to 
the applied pressure given in (a) for L = 100 with k = 0.02. The dotted line marked (r) indicates 
the long-time average of the period-averaged CO2 production rate r, defined by Eq. while 
the dotted line marked (i?c02)max marks the maximum average CO2 production rate for constant 
y = y2{k), shown in Fig. ^ Time is measured in units of MCSS. 

The transition line depends strongly on the selected values of yi and i/h, as can be seen 
by comparing Fig. IHl (a) and Fig. IHl (b). This behavior can be understood by looking at the 
dependence of the lifetimes on the pressure. In Fig. |3]we show a schematic, but qualitatively 
correct, diagram of the lifetimes Tp and r^vsy. The lifetimes increase rapidly as y approaches 
the coexistence point y2{k), but the dependences of and Tp on the distance to y2{k) are 
not equal (this is evident in the limiting case. A; — ^ 0, where ^ 00 independent of y), as 
the lifetimes depend on the direction of the decay. For the values of yi and yh selected in 
Fig. Ufa), the system takes longer to be decontaminated than to be poisoned, i.e., > Tp, 
suggesting that the transition line (corresponding to the region of high CO2 production) 
lies in the region > tp, as can be seen in Fig. El (a). When yi and yh are selected as in 
Fig.m^b), the system takes less time to be decontaminated than to be poisoned, i.e., < Tp. 
Then the transition line is expected to lie in the region td < tp, as seen in Fig. IHl (b). The 
location of the transition line in each case is schematically indicated in Fig. 01^c). 




1250 



1500 




7 





FIG. 3: Long-time average (r) of the period-averaged CO2 production rate r, shown as a density 
plot vs tfi and tp for (a) yi = 0.52, yh = 0.54, and (b) yi = 0.51, yh = 0.535. k = 0.01 and L = 100. 



In Fig. El it can be seen how the average CO2 production, (r) , depends on for two 
different values of tp and k. When k = 0.01, Fig. El (a) shows [consistent with Fig. El (a) 
and Fig. El (b)] that the system has two well-defined dynamic phases, one with (r) > and 
the other with (r) ^ 0. When k is increased to 0.04, the system changes continuously from 
one phase to the other, as seen in Fig. El (b) and Fig. El This behavior suggests that for 
low values of k, the system has a dynamic phase transition (DPT) between a high CO2 
productive dynamic phase and a nonproductive one. This is reminiscent of the behavior of 
ferromagnetic Ising or anisotropic XY or Heisenber g sp in systems driven by a periodically 



oscillating field |l8, 
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Universality and finite-size scaling are well-known tools to analyze critical phenomena. 
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FIG. 4: (a) and (b): Schematic representations of the decay times of the metastable states as 
functions of the CO pressure y. The asymmetry of the curves indicates that the hfetimes depend 
on the direction of approach to the coexistence point, y2{k). The ratio r^/rp depends on the values 
of yi and y^, which are different in (a) and (b). (a) corresponds to the situation shown in Fig.|21(a), 
and (b) corresponds to Fig.|21(b). (c) Schematic representation of the dependence of the transition 
hne between the phases with (r) > and (r) w 0, t*^{tp), on how yi and y^ are selected. 
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FIG. 5: Long-time average of the period-averaged rate of CO2 production, (r), shown vs for two 
values of tp and L = 100; (a) with yi = 0.52, yh = 0.535, and k = 0.01, and (b) with yi = 0.52, 
yh = 0.553, and k = 0.04. Only for k = 0.01 the system clearly presents two dynamic phases: one 
with (r) ~ and the other with (r) > 0. 
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FIG. 6: Surface plot of (r) vs td and tp for yi = 0.52, yh = 0.553, and k = 0.04. L = 100. 



Near a second-order equilibrium phase transition, the order-parameter correlation length 
diverges, leading to power-law singularities in terms of the finite system size L Q . Follow- 
ing the lead of prev ious applications to the DPT in the kinetic Ising model driven by an 
oscillating field j^, Q , we use the period-averaged CO2 production rate r as a dynamic 
order parameter and perform a finite-size scaling analysis of its fluctuations and mean value. 
We define a measure of the fluctuations in r in a L x L system in the standard way as 

X, = L^[{r^) - {r^] (3) 

and measure Xl as a function of td for a fixed value of tp and several values of L. Anal- 
ogous to the situation at a second-order equilibrium phase transition, the order-parameter 
fluctuations increase with the system size, such that the maximum value of Xl scales as 
^max ^ L"'/'^^ and the nth. moment of the order parameter at the transition scales as 
{r^)L ~ L""^^/*^^. (We use standard notation for the critical exponents: 7 for the fluctuation 
exponent, (3 for the order-parameter exponent, and u for the correlation- length exponent 

In Fig. [3 we show Xl vs td for four system sizes at A; = 0.01 and k = 0.04. The errors 
in Xl are obtained by standard error propagation analysis as ctx^ ~ '^Xl/^/ti — 1. Fig. [7| 
(a) shows that, for k = 0.01 and the four values of L used, Xl displays a clear peak, 
which increases in height with increasing L. This is in clear contrast with Fig. [7| (b) for 
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FIG. 7: The order-parameter fluctuation measure Xl, shown vs for tp = 230, for four system 
sizes, (a) k = 0.01 and (b) k = 0.04 w kc- The dotted hues are guides to the eye. For clarity the 
error bars are not included in the plot. The highest values of Xl have an error of approximately 
6 %. 



^ kc, which shows no such increase with L. 



k = 0.04 

In Fig.lSl^a) we plot ln(X™'^^) versus ln(L) for k = 0.01. A linear fit indicates a power-law 
divergence of the fluctuations with L, with exponent ^jv = 1.77 ±0.02. A different method 
to extract the power-law exponent, which has some advantage in eliminating the effects of a 
nonsingular background term (as in X/, = f + gL"'^'^ with / and g constants), is to consider 



In 



Xmax 
bL 
J^max 



/ln6 
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C(l/ln6) 



(4) 



with L fixed at a relatively small value (here, L = 80), and 6 > 1. For large L and b, 
the correction term is proportional to f/{glnb), so that the exponent can be estimated 
by plotting the left-hand-side of Eq. (0} vs l/ln6 and extrapolating to l/ln6 = 0, as in 
Fig-Efb). The resulting estimate is again 7/// = 1.77 ± 0.02. 

To obtain an estimate for the exponent ratio /?/z/, we used the scaling relation for the order 
parameter at the critical point, (r") ~ L~'^^^^/'^\ (again with the possibility of a nonscaling 
background, which can be quite significant for r) to plot the left-hand side of 



-In 



bL 



(r")j 



/In6 = n- + 0(l/ln6) 



(5) 
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FIG. 8: (a) Plot of ln(Xf vs ln(L) and (b) plot of ln(X^^^^/X]^^^)/ ln(6) vs l/ln(6), both for 
k = 0.01 and tp = 230, and including all four system sizes. Xf^^^ is the maximum value of X^, 
taken from Fig.[7|(a). The straight lines are the best linear fits to the data and in both cases give 

^max ^ j^^/u ^i^j^ ^ i ^jj ^ Q 

for L = 80 vs 1/ In 6, as shown in Fig. IHl For lack of a better estimate of the transition 
point, (r") were measured at the values of td corresponding to the maxima of Xl. A linear 
fit that takes into account all four system sizes gives (3/v = 0.14 ± 0.06 for n = 2 and 
(3/v = 0.10 ± 0.03 for n = 4. As a combined estimate, we take P/u = 0.12 ± 0.04, where 
the error bar includes some measure of our uncertainty about a nonscaling background and 
other finite-size effects. 

Combining the exponent estimates we find 



where d is the spatial dimension. This result agrees with hvperscaling and indicates 
that our finite-size scaling results are internally consistent [2^- Results very similar to the 
ones presented above were obtained selecting the period-averaged CO coverage as the order 
parameter, instead of r. 

Our estimated exponent ratios are very close to the analogous two-dimensional equi- 
librium Ising values {-f/u = 7/4 = 1.75 and p/u = 1/8 = 0.125 with /? = 1 Q), but 
they are also near those for two-dimensional random percolation (7/1^ = 43/24 ^ 1.79 and 



2{P/u) + (7/z/) = 2.01 ± 0.03 ^ 2 



d, 



(6) 
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FIG. 9: Plot of -ln((r")5i/(r")L)/ln(6) vs l/ln(6) for k = 0.01 and tp = 230, including all four 
system sizes. The straight lines are the best linear fits to the data, giving (3/u = 0.14 it 0.06 for 
n = 2 and P/u = 0.10 ± 0.03 for n = 4. 

(3/v = 5/48 ~ 0.104 with v = 4/3 ~ 1.33 One way to verify the universality class 

without calculating u directly (which would require much more accurate data than we have 
available), is to consider another universal quantity, such as the fixed-point value of the 



- 

fourth-order order-parameter cumulant ( "Binder cumulant" ) j30 1 



3((r-(r),)^)i' 

where denotes the average over the whole time series of r for an L x L system. This 
cumulant is shown vs td for different L in Fig. E3 The maximum possible value of is 2/3, 
which is reached in the dynamically ordered phase, provided that ergodicity is not broken 
or that (r) l is exactly known. This is not so in the present case, and so the cumulant 
is nonmonotonic, becoming negative deep in both the dynamically ordered (large td) and 
disordered (small t^) phases. With sufficiently accurate data, the curves representing ul for 
different L cross or touch at a common point, which represents an estimate for the critical 
value of the control variable (here, td) that is quite insensitive to corrections to scaling. The 
value of Ul at this fixed point, m*, is a universal quantity characteristic of the particular 
universality class. In the present case, the accuracy of our data is not sufficient to use the 
crossing point as an estimate for the critical value of td (we have instead used the maxima 
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FIG. 10: The fourth-order cumulant ul, shown vs td for all four system sizes, k = 0.01 and 
tp = 230. The horizontal lines correspond to ul = 2/3 (solid), tijsing ~ 0.610 (dashed), and 
''^Perc ^ 0.555 (dotted). See discussion in the text. The errors are approximated by standard error- 
propagation methods and are of the order of 3 % for the highest values of li^ and of the order of 
14 % for the lowest. 

of Xl). However, u* is seen to be in the vicinity of 0.6, consistent with the very accur ately 
known value for the two-dimensional Ising universality class, u^^- = 0.6106901(5) j3|. 
Even though it is not very accurately determined, the value of u* observed for the present 
system is unlikely to be as low as the value for random percolation, Up^^.^ ^ 0.555 m. 
Our total numerical finite-size scaling evidence thus points in the direction that the DPT in 
this system belongs to the two-dimensional equilibrium Ising universality class, together with 



other DPTs in far-from-equilibrium systems, such as the one o 



^served in the two-dimensional 
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25 



kinetic Ising model driven by an oscillating applied field 

An independent indication that the DPT in this system belongs to the equilibrium Ising 
universality class, is a symmetry argument due to Grinstein, Jayaprakash, and He jsTj, 
who argued that the equilibrium Ising universality class extends to nonequilibrium cellular 
automata that (i) have local dynamics, (ii) do not conserve the order parameter or other 
auxiliary fields, and (iii) respect the Ising up-down symmetry. This result was later extended 
by Bassler and Schmittmann jssi, using renormalization-group arguments, to systems that 
obey conditions (i) and (ii) above, but violate (iii), such as some driven lattice gases. The 
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ZGB-k model satisfies these requirements, and it should thus belong to the equilibrium Ising 
class on symmetry grounds. 



IV. CONCLUSIONS 

In this paper we have studied the dynamic response of the ZGB model with CO desorption 
(the ZGB-k model jl^) to periodic variations of the relative CO pressure y, around the 
coexistence value y2{k) that separates the low and high CO-coverage phases. There is an 
asymmetry of the lifetimes of the model: its decontaminating time generally differs from 
the contaminating time Tp. We exploited this fact by selecting a square-wave periodic CO 
pressure that stays for a time td in the high-production region and for a time tp in the 
low-production one. We found that td and tp can be tuned to significantly enhance the time- 
averaged catalytic activity of the system beyond its maximum value under constant-pressure 
conditions - a result we believe should be of applied significance. 

We also found strong indications that, for sufficiently low values of the desorption rate, 
this driven nonequilibrium system undergoes a dynamic phase transition between a dynamic 
phase of high CO2 production, (r) > 0, and a nonproductive phase (r) ^ 0. As the order 
parameter for this nonequilibrium phase transition we used the period-averaged rate of CO2 
production, r. Our study shows that the distinction between these phases disappears for 
a high enough desorption rate. Applying finite-size scaling techniques in a similar fashion 
to what is commonly done to study equilibrium second-order phase transitions, we found 
that, for small values of the CO desorption rate fc, the fluctuations of the order parameter 
diverge as a power law with the system size, 

j5^max _ j^^/v ^-^^^i expoucut 7/// = 1.77 ±0.02, 
while moments of the order parameter at the transition point decay as (r") i ~ ij-'^P/^ with 
P/v = 0.12±0.04, and the fourth-order order-parameter cumulant has a fixed-point value 
of u* ~ 0.6. These values are close to those of the two-dimensional Ising universality class, 
and together with general symmetry arguments js? , they represent reasonable evidence 
that this far-from equilibrium phase transition belongs to the same universality class as the 
equilibrium Ising model. A a higher level of confidence about the universality class would 
require simulations at least an order of magnitude more extensive than the ones presented 
here - a task beyond the scope of the present paper. 

Finally, we note that the enhancement of the CO2 production and the continuous DPT 
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are probably closely related phenomena since the critical cluster associated with the phase 
transition is likely to provide more empty sites available for O2 adsorption near adsorbed 
CO molecules, than the sharp interfaces expected near the first-order coexistence line seen 
under constant-|/ conditions. 
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